#ifndef FEHLBERG_HPP
#define FEHLBERG_HPP

#include "EmbRK.hpp"

class Fehlberg : public EmbRK {
public:
  Fehlberg(int dim, double dt, double t0, double t1) : EmbRK(dim, dt, t0, t1) {
    m_s = 6;
    m_ph = 5;
    m_pl = 4;

    if (!initializeSolver())
      std::cerr << "Error in ClassicalRK initialization." << std::endl;
  }

  bool initializeSolver() {
    m_A.resize(m_s, m_s);
    m_c.resize(m_s);
    m_bl.resize(m_s);
    m_bh.resize(m_s);
    
    // 初始化系数
    m_A.setZero();
    m_c.setZero();
    m_bl.setZero();
    m_bh.setZero();
    
    // 填充 c 向量
    m_c << 0.0, 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0;
    
    // 填充 A 矩阵 (下三角部分)
    m_A(1, 0) = 1.0/4.0;
    
    m_A(2, 0) = 3.0/32.0;
    m_A(2, 1) = 9.0/32.0;
    
    m_A(3, 0) = 1932.0/2197.0;
    m_A(3, 1) = -7200.0/2197.0;
    m_A(3, 2) = 7296.0/2197.0;
    
    m_A(4, 0) = 439.0/216.0;
    m_A(4, 1) = -8.0;
    m_A(4, 2) = 3680.0/513.0;
    m_A(4, 3) = -845.0/4104.0;
    
    m_A(5, 0) = -8.0/27.0;
    m_A(5, 1) = 2.0;
    m_A(5, 2) = -3544.0/2565.0;
    m_A(5, 3) = 1859.0/4104.0;
    m_A(5, 4) = -11.0/40.0;
    
    // 填充低阶 (4阶) 系数
    m_bl << 25.0/216.0, 0.0, 1408.0/2565.0, 2197.0/4104.0, -1.0/5.0, 0.0;
    
    // 填充高阶 (5阶) 系数
    m_bh << 16.0/135.0, 0.0, 6656.0/12825.0, 28561.0/56430.0, -9.0/50.0, 2.0/55.0;

    // m_b = m_bl;
    
    return true;
  }

  // Eigen::VectorXd m_bl;
  // Eigen::VectorXd m_bh;
  // int m_pl;
  // int m_ph;
};

#endif